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Abstract 

For the diffraction of an incident plane electromagnetic wave by a slotted metallic film, the 
previous analytical calculation for a single slot [Technical Phys. 50, 1076 (2005)] is generalized 
into a model for an arbitrary linear array of slots with variable slot width, slot separation and slot 
dielectric material. The advantage as well as the effectiveness of the generalized model presented 
in this paper are best described by enabling calculation of a continuous spatial distribution of 
an electromagnetic field by inverting a small discrete coefficient matrix spanned by both the slot 
index and the slot-eigenmode index for a set of linear equations. In comparisons with well-known 
plane-wave and finite-difference time-domain methods, inverting a large matrix, in wave number 
space for the former case and in real space at each time step for the latter case, can be avoided 
to greatly speed up numerical calculations. In addition, based on a partial-domain method, the 
formalism presented here can be employed to treat a composite surface (e.g., a slotted metal film 
with different dielectric materials in the slots), while the analytical Green's function approach [J. 
Opt. A: Pure Appl. Opt. 8, S191 (2006)] becomes intractable in this case. Some numerical results 
are presented as a demonstration of this new analytical model. 



I. INTRODUCTION 



For a periodically-patterned surface on a metal film, the Maxwell equations can be 
solved by using a plane-wave (PW) method,'^ or a finite-difference time-domain (FDTD) 
method,'^ or a Green's function (GF) method.'^ For the PW method, a large matrix in 
wave number space with respect to different reciprocal lattice vectors needs to be inverted 
in order to obtain the spatial distribution of the total electromagnetic (EM) field. For the 
FDTD method, the time evolution into a steady state EM field is obtained through in- 
verting a large matrix in real space at each time step with respect to all grid points in a 
region containing a unit-cell and a surrounding space not far from the surface. For the GF 
method, all the values of the EM field on a smooth surface must be calculated by inverting 
a large matrix in real space with respect to all grid points on the boundaries to get the 
spatial distribution of the total EM field. Both the PW and GF methods have an advantage 
for calculating a far-field distribution, while the FDTD method is usually limited to the 
calculation of the near-field distribution. Moreover, the GF method can also be used for a 
non-periodic smooth surface profile function. 

For a non-smooth surface, on the other hand, such as a one-dimensional periodic slot 
array, the GF method cannot be used. In this case, however, the partial-domain method'^ 
used in the GF theory can still be employed in combination with the slot-eigenmode ex- 
pansionl^ as well as with the Fourier expansion'^El methods. Based on this technique, an 
analytic solution for the diffraction of a plane EM wave by a single slot has been obtained.'^ 
This analytical single-slot approach has been generalized in this paper to deal with an arbi- 
trary slot array. If a single slot is replaced by double slots, each surface-plasmon-polariton 
branch'HHIS ^ffl[\ Iqq split into two with a minigap controlled by an EM coupling between the 
two slots. In the case of double slots with different widths and filled with different dielec- 
tric materials, the circulation and weaving of light was found as a result of the excitation of 
a phase resonance.'^ 

In this paper, a two-dimensional (2D) theoretical model for the diffraction of a plane EM 
wave by a perfectly electrically conducting (PEC) slotted film is derived with variable slot 
width, slot separation and slot dielectric material. The solutions of the 2D Helmholtz equa- 
tion in different partial domains are matched to each other by proper boundary conditions 
for both PEC and slot dielectric medium. This method has applications in modeling surface 
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plasmon polaritrons or surface waves, '^optical beam steering, and metamaterial design. 
The use of a PEC material precludes modeling of dissipation but allows for simple modal 
expressions for the slot fields. '^'^ When the optical depth of a metallic film is small and 
energy dissipation is not a consideration (off-resonance), the use of a PEC partial domain 
model is expected to be a good approximation. 

This paper is organized as follows. In Sec.[TT| we generalize the single-slot model into 



one for an arbitrary array of slots with both p and s polarizations. In Sec. Ill, numerical 
results as a demonstration are presented and compared for the transmitted EM wave with 
p polarization through a slot array with different slot numbers and incident angles. A brief 



summary of the generalized model is given in Sec. IV 



II. THEORY 

In this paper, the incident EM field is assumed to be a plane wave with a simple time 
dependence proportional to e"*"^*, where u is the angular frequency of the incident field. 
For the s-polarized case, we write E = (0, Ey, 0) = (0, m, 0) for the transverse electric 
field and H = [—i/{ujjio)] ^ x E. For the p-polarized case, on the other hand, we write 
H = (0, Hy, 0) = (0, M, 0) and E = [z/(a;eoe)] ^ x H. Here, eo and /io are the vacuum 
permittivity and permeability, c = l/i/eo /io is the speed of light in vacuum, and e is the 
relative dielectric constant of a non-magnetic dielectric medium. 

The master equation for determining the scalar function u introduced above is the two- 
dimensional Helmholtz equation in the dielectric medium in the xz— plane: 

_ + _+i=„ = 0, (1) 

where k = {u/c)^/e is the wave number of the field. For a perfect electric conductor (PEC) 
with an extremely-large conductivity, the boundary condition for s polarization associated 
with Eq. ([T| is just u = along all PEC boundaries. In addition, for p polarization, the 
boundary condition of Eq. ^ is du/dn = along all PEC boundaries where n is in the 
surface normal direction. 
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A. s— Polarization 



For the system shown in Fig.[T| Region I, below the slotted layer {x < —d), has the 
boundary condition 



u 



'd+0 ' 



for upper PEC 
for middle slot 
for lower PEC 



Zj + ij < Z < Zj+i 



Zj\ < Ij 



(2) 



: 

u\ 

, for lower PEC zj^i + iij^i < z < zj 

where j = 1, 2, ■ ■ ■ , is the slot index, Zj is the coordinate of the jth slot center, 2ij is the 
width of the jth slot, and 2d is the depth of all slots or screen thickness. For j = N, we take 
^AT+i = oo and In+i = in Eq. Similarly, zq = — oo and £o = in Eq. ^ when j = 1. 
Moreover, from the continuity of the derivative of u, we have another boundary condition 



du 
dx 



x=-d-0 



du 
dx 



for all middle slots 



\Z-Zj\ < 



(3) 



In a similar way, we are able to get the boundary conditions for Region III above the slotted 
layer (x > d) from Eqs. ([2]) and (jsj) by setting x = (i ± 0. Finally, the field must be zero 
along the interior slot walls in Region II (|a;| < d), that is. 



OC J Z — Zj —I ^ j J 







at all slot walls Ixl < d . 



(4) 



In Region I (below the slot layer), the solution of the Helmholtz equation in Eq. ([T]), 
including the incident field (incident from below), can be written as 

U^^\x z) = Q^^o{x+d) ^il3oz _ ^-iko(x+d) ^iPoz 
oo 

+ jdP [AsiP) cosiP z)+tAa{l3) sin(/3^)] e-*'=i('^)("+'^) , (5) 



where the subscripts s and a represent the symmetric and antisymmetric contributions, 
separately, ko = [uni^/c) cos 6*0, /3o = (i^'^l/c) sin 6*0, /? is real, /ci(/3) = ^J [uni^/cY — 
can be either real or complex with lm[A;i(/3)] > 0, til = ^/e[ is the refractive index for the 
dielectric medium below the slot layer, and is the incident angle. In Region III (above 
the slot layer), on the other hand, the solution of the Helmholtz equation takes the form of 
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df3 [Bs{f3) cos{/3z) +iBa{P) sm{/3z)] e 



ik2{l3)(x-d) 



(6) 



where = ^/eR is the refractive index for the dielectric medium above the slots and 
^2(/3) = a/ (ooTiji/c)^ — can also be real or complex with lm[fc2(/3)] > 0. 

In Region II (the slot layer), the solution of the Helmholtz equation along with the 
boundary condition in Eq. Q is found to be 



\z - Zj 



1)5: {[< 
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^(d+x) _|_ y ^iain{d-x) 



COs[^l^iz - Zj)] 



(7) 



where Q(^z) is the unit step function, ri = 1, 2, ■ ■ ■ is an integer for the eigenmode index. 



{n - 1/2) and 



in-n jlj^ are for the symmetric and antisymmetric slot 



modes, respectively, a- 



sn, an 



[urij/cy — (^sn, an)^ cau bc either real or complex with 



Re[cr^„^„] > 0, and rij = ^/e] is the refractive index for the dielectric medium inside the jth 
slot. 

By using the boundary conditions in Eq. ^ for u at x = ±d and the orthogonality of the 
Fourier expansions in Eqs. ([s]) and ([6]), the coefficients for both symmetric and antisymmetric 
contributions can be expressed by the slot eigenmodes as follows: 



Asm 



TT 



dzu^^^\x = —d, z) cos(/3z) 



5: El 



2iai„d 



- I f<„ + Kne''"'-') Qim sin(/3^,) 



(8) 



oo 

— / dzu^^^\x = —d, z) sin(/3z) 
J 



n K. J 



(9) 
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IT 



dzu^^^\x = d, z) cos{Pz) 



EE 

n L j 



(X „^ €■ 



(10) 



dzu^^^\x = d, z) sin(/3z) 



EE 



(11) 



In Eqs. we have defined the notations 



dz cos(/32;) cos(^^„2;) 



sinc[(/3-6J^,]+sinc[(/3 + ej^.] 



(12) 



<5L(/3) = Y- j sm{l3z) sin(^;^„2;) 

= sinc[(/3 - a) - sinc[(/3 + ^L) i,] , (13) 

where sinc(x) = sin(x)/x. 

The slot eigenmode expansion coefficients a^„, and in Eq. ([7| are determined 

from the derivative continuities in Eq. ([3]) for n at x = ± d and for each slot. Therefore, we 
get for J = 1, 2, ■ ■ ■ , iV 
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dz 
dz 



du^ix, z) 



dx 



du^^{x, z) 



x=-d-0 



dx 



(14) 



x=-d+0 



dz 
dz 



COs[^i^{z - Zj)] 
COs[^i^{z - Zj)] 



du^ {x, z) 



dx 



du^^{x, z) 



x=d+0 



dx 



x=d-0 



(15) 



Using the orthogonality of the slot eigenmodes, from Eqs. (14) and (15) we arrive at the 
following set of equations: 



2koe'^°^^QiM) - / dp [A,{P) cos{f3zj)+tA,{P) sm{f3zj)] QiMHP) 



(16) 



- 2fcoe*^°^^- Qi„(/3o) + / dp [iA,{P) sin(/3^,) + A,(/3) cos(/3z,)] Q{^{P) h{P) 



~<yL (aL - V.^ e^"^-") , (17) 



an \ an an 



dp cos(/3z,-) + iB^{P) sin(/3^,)] Qi„(/3) k^{P) = {ai^ e^-^-^^ - fe^J , (18) 



dp [iBsiP) siniPzj) + BaiP) cosiPz,)] QUP) hiP) 



"an I 1 



(19) 
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where j = 1, 2, ■ ■ ■ , A^. 



Inserting Eqs. (^8|-(11) into Eqs. (16)-(19) gives rise to a set of AN linear equations with 



respect to al„, and b{^. That is, for j = 1, 2, 



2ia^ ,d 



n',j' 



E 



^an' ' an' ^ 



TX 







2ia^ ,d\ 



(20) 



Ce"-^ 



TT 



an' 



dPQUP) Qilm km sin[/3(z, - zy)] 



dPQUP)QL'{P) km cos[(3iz,-z,,)] 



n',j' 



7r 



n',j' 



d^QLW) QUP) k2W) cos[/3(z,, - zj)] 



ioo 



QU^) Q^m k2W) sm[^{zy - z,)] 



n',j' 



(21) 



(22) 



{oo 


(oo 
n 



d^QirMQiUP)k2W) sm[Piz,-z,,)] 



dp QiM Qi'n'W) k2W) C0S[/3(Z, - Zy)] 



(23) 



If we truncate the number of slot eigenmodes in the expansion, i.e., n = 1, 2, ■ ■ ■ , M, this 
constitutes a 4iVM-order inhomogeneous hnear-matrix equation, which can be exactly solved 
by inverting its coefficient matrix. Here the symmetric and antisymmetric eigenmodes in 
different slots with Zj 7^ Zji are coupled to each other, which is different from the single-slot 
result.'^ 



B. p— Polarization 

Following an approach that parallels the approach we performed above to deal with the 
boundary conditions for s polarization, we require for p polarization that the derivative 
condition is interchanged with the null conditions, i.e., m = becomes du/dn = along the 
PEC surfaces. In addition, the {l/e)du/dx continuity condition across the slot interfaces 
is used to determine the Fourier coefficients outside Region II instead of the u continuity 
condition. For Region I (x < —d), we have the boundary condition 



1 du 
eL dx 



x=-d-0 



1 du 

e, dx 



for upper PEC zj + ij < z < Zj^i — £j+i 
for middle slot 



x=-d+0 



\Z Zj\ < Ij 



(24) 



for lower PEC Zj-i 



< z < z 



'3 



For Region III (a; > d), the boundary condition is 



1 du 



x=d+0 



1 du 
e, dx 



x=d-0 



for upper PEC zj + ij < z < zj^i — ij^i 
for middle slot \z — Zj\ < ij 

for lower PEC zj^i + ij^i < z < zj — ij 



The continuity of u is required along the interfaces of each slot: 



(25) 



^U=-d-0 '^\x=-d+0 ' 



'^\x=d~0 ~ ''^\x=d+0 ' 



for all middle slots \z — Zj\ < ij 
for all middle slots \z — Zj\ < ij 



(26) 



where j = 1, 2, ■ ■ ■ , A^. Finally, the field normal derivative must be zero in Region II along 
the interior slot walls (|x| < d), that is, 

du{x, z = Zj ± ij) 



dz 



, at all slot walls \x\ < d 



(27) 
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In Region I, the solution of the Helmholtz equation in Eq. ([T]) in this case, including the 
incident field, can be written as 



dl3 



[As{(3) cos(/3z) + iAa{(3) sin(/3z)] e-*'=i(^)(^'+'^) . (28) 



In Region III, the solution of the Helmholtz equation takes the form of 



dl3 

km 



[Bs{(3) cos{(3z) + iBaiP) sm{(3z)] e*'=2(/3)(x-d) _ ^29) 



In Region II, the solution of the Helmholtz equation along with the boundary condition in 



Eq. (27) is found to be 



mW(x, z) 



11. 

sn 



_|_ I j Li ^i<^irzi.d+x) _ ^i(Ti„(d-x) 



COs[^i^{z - Zj)] 

(30) 



where = (n/ij) {n — 1) and = {ji /dj) {n — 1/2) are for symmetric and antisymmetric 
slot eigenmodes, respectively. 



By using the derivative boundary conditions in Eq. (24) and (25) at a; = ± ci as well as the 



orthogonality of the Fourier expansions in Eqs. (28) and (29), the coefficients As(/3), v4a(/3), 
Bs{f3) and Ba{f3) can be expressed as 



AsiP) = - dz- ^— 

— CO 



cos{f3z) 



x=—d 



E 1 E 7 [{<n + Kne''''-') QinW) COs(/3z,) 



QUP) sin(/3^,)]} 



(31) 
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lAaiP) = - / dz — 

— oo 



sm{l3z) 



x=—d 



5: 



Qi„(/3) sin(/3z,) 

/I L \ / 



(32) 



1 f ^ 1 

— oo 



cos(/32;) 



x=d 



n L j 



— I \ a 



(33) 



iBam = - dz- ^— 

— oo 



sm{(3z) 



x=d 



(34) 



The slot eigenmode expansion coefficients a^„, and 6^^ in Eq. (30) are determined 



from the continuities of u in Eq. (26) a.t x = ±d for individual slot. As a result, we get for 
J = 1, 2, ■ ■ ■ , iV 



Zj +£j 



dz 
dz 



cos[^f„(2; - Zj)\ 



^3 *-j 



u^(a; = — (i, 2;) 



cos[Cf„(2; - 2;^)] 



u^^ix 



-d, z) , 



(35) 
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dz 



u {x = d, z) 



dz 



sin[^L(^ - Zj)] 



u^\x = d, z) . (36) 
As for the s— polarization case, using the orthogonahty of the slot eigenmodes, as well as 



using Eqs. (31)-(34), gives rise to the following set of AN linear equations with respect to 



^sn^ ^sni ^an ^^'^ ^an 



TT 



d/3QU/3)Qi^,iP) 



1 



km 

1 



km 



co?,[[5{zji — Zj)] 
sin[/3(zj/ - Zj)] 



I \ ^ ^j,j'^n,n'Xn'(^j' / f ti' 
+ 2^ j' -Kn'-Kn'^ 



n',j' 



2ia-' ,d\ 

sn' I 



ko 



Jl3oZj r~)j 



QM) 



(37) 



{oo 




A;i(/3) 



dPQinW)Qi'nm 1-7^ ^OsiPiz, - Z,>)] 



km 



ko 



QU(3o) 



(38) 



TT 



TT 



^ dp QiM qUP) MPizy - z,)] 



k2W) 



0"; 



(39) 
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oo 



n',j' ( ^ 

{oo 


_ 5^^M:^(,/^,e-<>-6L) = 0. (40) 

Here, Xn = 2 for n = 1 and Xn = 1 for ^ 7^ 1- Similarly, if we limit the number of slot 
eigenmodes in the expansion to M, this constitutes a 4A^M— order inhomogeneous linear- 
matrix equation which can be solved exactly by matrix inversion. Again, the symmetric and 
antisymmetric eigenmodes in different slots with Zj 7^ Zj' are coupled to each other. 



III. NUMERICAL RESULTS 



To demonstrate these results, we perform numerical calculations here for the 
polarization case which allows for the excitation of a surface-plasmon-polariton mode. In 
our numerical calculations, we have taken: Aq = 27tc/uj = 1.0 fim (for a large transmission), 
£j = Ao/10 = 0.1 /im (deep sub-wavelength regime), d = 3.0 £j = 0.3/im, riL = ""-r = = 1, 
and the center-to-center slot separation is assumed to be 0.98 fim. Other parameters, such 
as ^0 and A^, will be directly given in the figure captions. 

Figure [2] displays the three-dimensional plot of the transmitted EM wave \u{x, z)p = 
\Hy{x, z)p in Region III for p polarization and normal incidence with N = 2 and ^0 = 0°. 
Here, the strong transmitted near-field close to the exits of the two slots in the z direction 
decay very fast as it propagates away from the slot array in the x direction. At the same 
time, a local maximum is gradually built up in the center region between the two slots due 
to strong interference of these two slot-exit waves and it becomes weaker away the slot layer 
in the x direction. When the number of slots is increased from two to five in Fig.|3j the 
strong transmitted near-fields close to the five slot exits interfere with each other strongly 
in the z direction and constitute a complex pattern as they propagate away from the slot 
array in the x direction. There are four local maxima that develop eventually at the four 
centers between two adjacent slots as the strong slot-exit near- fields are partially suppressed 
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with distance away from the array in the x direction. For normal incidence, although the 
symmetric and antisymmetric modes of different slots are coupled to each other, the mirror 
symmetry with respect to the center of the system (at z = 0) still remains. However, this 
mirror symmetry is broken for a nonzero incident angle, as demonstrated in Fig.|4]by the 
two-dimensional contour plot of the transmitted EM wave \u{x, z)p = \Hy{x, z)\'^ in Region 
III for p polarization with = 60°. In this case, the transmitted EM wave becomes weaker 
in general. At the same time, the field strength increases with z from the lowest slot to the 
upper most slot, which is accompanied by an increasingly extended spatial distribution of 
the fields between two adjacent slots as one moves up in the z direction. 

The advantage of the current approach is in casting a continuous spatial distribution 
problem into a small discrete coefficient matrix spanned by the slot index j and the eigen- 
mode index n. In this paper, we only display the spatial distribution of a transmitted EM 
wave in the deep sub-wavelength regime [ij <^ Aq). However, the accuracy of our numer- 
ical calculations becomes worse as the slot width 2£j becomes comparable to the incident 
wavelength Aq. In this case, a very large number of slot eigenmodes are required to achieve 
a good accuracy in a numerical calculation. 

IV. CONCLUSIONS 

In conclusion, the previous analytical calculation for the diffraction of an incident plane 
electromagnetic wave by a single slot on a metal film has been generalized to treating an 
arbitrary linear array of slots with variable slot width, slot separation and slot dielectric ma- 
terial. In comparisons with the plane-wave and finite-difference time-domain methods, the 
calculation speed has been increased greatly by transforming the calculation of a continuous 
spatial distribution of an electromagnetic field to the inversion of a small discrete coefficient 
matrix spanned by the slot and eigenmode indexes. Moreover, based on a similar partial- 
domain method, the linearly slotted metal film with different dielectric material in each slot, 
which cannot be solved by using the analytical Green's function approach, has been studied 
by our analytical approach. Some numerical results that demonstrate our generalized model 
have also been presented. 
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FIG. 1: (Color online) Illustration of a one-dimensional slot array in the xz— plane, where j = 
1,2, • • • , is the slot index, zj is the center of jth slot, 2ij is the jth slot width, the slot depth is 
2d, and is the incident angle. The regions below and above the slot array are defined as Region 
I and Region III, respectively. The slot-array region is designated Region II in the text. 
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FIG. 2: (Color online) Three-dimensional (3D) plot of transmitted EM wave \u{x, z)p = 
\Hy{x, z)p in Region III for p polarization with N = 2 and = 0° for normal incidence. 
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FIG. 3: (Color online) 3D plot of transmitted EM wave \u{x, z)p = \Hy{x, z)p in Region III for 
p polarization with N = 5 and = 0° for normal incidence. 
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